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Abstract 

The Schrodinger equation for two and tree-body problems is solved for scattering states in a hybrid repre- 
sentation where solutions are expanded in the eigenstates of the harmonic oscillator in the interaction region 
and on a finite difference grid in the near- and far-field. The two representations are coupled through a 
high-order asymptotic formula that takes into account the function values and the third derivative in the 
classical turning points. For various examples the convergence is analyzed for various physics problems that 
use an expansion in a large number of oscillator states. The results show significant improvement over the 
JM-ECS method [Bidasyuk et al, Phys. Rev. C 82, 064603 (2010)]. 
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1. Introduction 

The accurate prediction of the breakup of a many-particle system into multiple fragments is one of the 
most challenging problems in quantum mechanics. Not only the relative motion of the particles needs to be 
modeled, but also the internal structure of the target and the products need to be described accurately. This 
leads in many cases, to a high-dimensional Schrodinger equation posed on a huge domain. For example, in 
breakup or collisions of nuclear clusters the cross section depends on a delicate interplay of the forces that 
hold the clusters together and the forces between the clusters. 

The quantum state that describes the internal structure of a bound many particle system is often represented 
as linear combination of eigenstates of the quantum mechanical oscillator. These states form an L^-basis 
which reduces the problem to finding the correct expansion coefficients. Examples are the correlated many 
electron state of a quantum dot [2 El E] and the many nucleon state in nuclear physics [H [5] . The various 
applications of the oscillator representation are discussed in [B] . Such a representation is very efficient for 
tightly bound ground state or lowest excited states as any smooth potential well is close to parabolic shape 
near it's bottom. As a result, if the oscillator parameter is optimized to match this local parabolic potential, 
the lowest states of a system can be efficiently represented with only a few oscillator functions [7] . 

However, the representation is inefficient to describe scattering or break-up processes. Scattering states are 
not square integrable and many oscillator states are required to represent the interaction and asymptotic 
region. Furthermore, many potential matrix elements need to be calculated and this results in a dense linear 
system that has a complexity of to arrive at the solution, where N is the number of oscillator states used 
in the representation, omitting the cost of calculating the potential matrix elements. 
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The J-matrix method offers a way of calculating cross sections and other scattering observables in the 
oscillator representation. It was proposed by Heller and Yamani in the seventies [3 [9] and mainly applied 
to atomic problems. The method exploits the tridiagonal structure of the kinetic energy operator. The 
J-matrix has been under constant development since its inception and a review of the recent developments 
can be found in [TU]. For problems with Coulomb interactions the Coulomb- Sturmian basis is preferred 
over the oscillator representation. Recently the Coulomb-Sturmians have been used to describe multiphoton 
single and double- ionization [11 and electron impact ionization |12j . 

At the same time, and mostly independent, the Algebraic Resonating Group Method was developed for 
nuclear scattering problems [H [T31 [TU [TSl US] . This method exploits the same principles as the J-matrix 
method for description of nuclear cluster systems, where the oscillator representation is efficiently used to 
describe the internal structure of clusters. If the same representation is used for intercluster degrees of 
freedom, then the nucleon symmetrization rules become straightforward. 

One shortcoming of the methods based on an basis is that the asymptotic solutions need to be known 
explicitly before a system for the wave function in the inner region and the scattering observables can 
be written down. This is a serious limitation since it is hard to find the asymptotic wave function for 
breakup reactions with multiple fragments. The Complex Scaling method, which scales the full domain into 
the complex plane, can be easily implemented in the oscillator representation by taking a complex valued 
oscillator strength. This has been used to calculate energy spectrum or extract the resonances [T7], but the 
calculation of cross sections and other scattering observables may be quite difficult. 

In recent decades significant progress has been made in the numerical solution of scattering processes de- 
scribed by the Helmholtz equation. In contrast to many particle systems, where the potential V is often 
non-local, the wave number k{x) in the Helmholtz equation depends only on the local material parameters 
such as the speed of sound in acoustics or the electric permittivity and magnetic permeability for electro- 
magnetic scattering. For these problems grid based representations such as finite difference, finite element 
[T5J [in] , or Discontinuous Galerkin pOl [5T] are preferred since they lead to sparse matrices that can be easily 
solved by preconditioned Krylov subspace methods [221 [13] . 

Another technique that has found widespread application is the use of absorbing boundary conditions. These 
boundary conditions allow a scattering calculation without prior knowledge about the asymptotic wave form. 
Exterior Complex Scaling (ECS) and Complex Scaling (CS) are widely used in atomic and molecular physics 
[231 [23 [23 [23 ■ Perfectly matched layers (PML) are used for electromagnetic and acoustic scattering [25] , 
which can also be interpreted as a complex stretching transformation |29] . There are many other excellent 
absorbing boundary conditions [30l [311 [32] . 

In [33] the JM-ECS method was introduced that combines the J-Matrix method with a grid based ECS. The 
method describes the scattering solution in the interior region with an oscillator representation and in the 
exterior region with finite differences. The two representations are matched through a low order asymptotic 
formula with an error that scales as 0[N^^/'^)^ where N is the size of the oscillator basis describing the 
inner region. Once the grid and oscillator representation arc matched, it is easy to introduce an absorbing 
boundary layer since the grid representation can be easily extended with an ECS absorbing layer or any 
other absorbing boundary condition. 

The resulting method was illustrated for one- and two-dimensional model problems representative for real 
scattering problems with local interactions. Furthermore, the representation was used for nuclear p-shell 
scattering. However, the accuracy of the calculations was unsatisfactory due the low-order matching condi- 
tion. While the grid representation on the exterior has an accuracy of 0{N~^), the matching is only accurate 
to order 0(N~^/^). 

The main contribution of the current paper is to increase the accuracy of the asymptotic formula that allows 
a better matching of the grid and oscillator representations. The better asymptotic formula takes function 
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values and its third derivative into account. It can bring the matching error down to the level of the accuracy 
of the grid representation. 

A higher-order asymptotic approximation of the oscillator representation was already discussed by S. Igashov 
in |34| . However, the formula was not used to increase the accuracy of scattering calculations. 

This paper is outlined as followed. In Section [2] we shortly describe the process of scattering calculations 
in the oscillator representation. In Section |3] we derive a higher-order asymptotic formula that takes into 
account the behavior of the function in the classical turning points in coordinate and Fourier space. In 
Section 4 and 5 we use this asymptotic formula to solve scattering problems. 



2. Review of scattering calculations in the oscillator representation 

In this section we discuss the most important properties of the oscillator representation and how they can be 
used to perform scattering calculations with the J-matrix method. We also recall the working of the hybrid 
J-matrix and ECS method proposed in |33j . 



2.1. The radial scattering equation 

The aim is to solve the Schrodinger equation in atomic units (m ^ 1, h — 1) that describes a scattering 
process of two particles for an energy i? G M. This equation written in relative coordinates is 



-^-A + Vir)-E 



*(r) = F(r), Vr e 



(1) 



where F{r) is the function describing the initial state or the source term and V{r) is the potential. The 
coordinate r can be written in spherical coordinates {p,6,(j}). In case of spherically symmetric potential 
(T^(r) = V{p))^ Eq. ([I]) can be reduced to a one-dimensional radial equation using partial wave decomposition 



l.m 



p 



F{r) = F{p,e,, 



p 



where Yi_rn{0,4') is a spherical harmonic, I = 0,1,2,... is the orbital angular momentum of the relative 
motion, m is the projection of this angular momentum. The resulting reduced radial equation becomes 

1 



2p2 



Mp) = Mp)- 



(2) 



For the problem we are interested in there is a range a > such that for all p > a both V{p) and x[p) are 



zero. 



We solve the Eq. ^ for E > Q and extract from the solution ?/'/ scattering observables such as the cross 
sections or phase shift. In order to solve the Eq. ([2| we represent the solution as 



(3) 



where tfn,i{p) are orthogonal functions, in particular we will use reduced oscillator functions, whose 
properties will be explained in the next section. After projection of ([3| on (pn,i, Eq. ([2| results in an infinite 
linear system 



E]c„ 



bk, 



(4) 
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where denotes the elements of the kinetic energy matrix and Vkn the elements of the potential energy 



matrix. They are the integrals 







and bk.i — Vk.i{p)xi{p) dp. As the considered problem is effectively one-dimensional we will use x instead 
of p to denote radial relative coordinate in all following sections devoted to two-body problem. 

2.2. The oscillator representation 

Before we explain the strategy to solve Eq. Q , we repeat the main properties of the oscillator representation 
and the function (pn^i that will be used in the expansion. 

The reduced radial equation analogous to Q for the quantum harmonic oscillator with an oscillator strength 
u) is 

r 1 ^2 ^ U] \ 1 _ _1 

(Pn,l{x) = En,liPn,l{x) (5) 



1 (P 11(1 + 1) 1 o n 



with x € [0, oo[ a radial coordinate. The boundary conditions are ipn,iiO) — and lim^^-s-oo 'fin.i{x) = 0. The 
eigenvalues are 

En,i = (in + ? + I ) (6) 



2 

and eigenstates 

„2 



62 



(7) 



where L^rt^^"^ are Laguerre polynomials. The normalization is iV„.j = ^/2n\/T(n + I + 3/2), where n € 
{0, 1, . . . , } and oscillator length is defined as 6 = -y/l/cj. 

The classical turning point associated with each state is Rn.i = b\/An + 2/ + 3, defined as the point where 
the potential energy equals the total energy of the system. 

The functions ([t]) form a complete orthonormal basis and any wave function il)i{x) that behaves as in 
X = Q can be represented using the infinite sum: 

/■oo 

= ^ c„,iV9„,;(a;), where c„,i = / ipn,i{x)^Ji{x)dx. (8) 

n=0 -^0 

In the following section we will use that the radial oscillator equation, ([5]), can be rewritten in terms of b, 
the oscillator length, and Rn,i, the classical turning point as 



^nA^)=0. (9) 



2(ix2 x^ 264 264 
An important property that forms the basis of the results in this paper is that the n-th oscillator state has 



n oscillations between the origin and the classical turning point i?„_; = y/An + 21 + 3. Beyond this turning 
point the function is exponentially decaying without additional oscillations. This means that as n increases, 
the frequency of the oscillation between the origin and Rnj grows proportional to y/n. This property will 
be used to derive the asymptotic formula in Section [3] 

The Bessel transform Fi{k) of a function Fi{x) is defined as 

Fiik) I Fi{x)n{kx)dx. (10) 
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Figure 1; The reduced oscillator state ip„j{x) with a classical turning point Rn,i for I = and two values of n: n = 10 (left) and 
n = 30 (right). The function oscillates n times on the interval [0, Since R„ ; only grows with 0{y/n) the state becomes 

rapidly oscillating as n grows. 

where ji{kx) is a Riccati-Bessel function, the regular solution of the radial Schrodinger equation without 
potential, ([2]). It is connected with the ordinary Bessel function in a following way ji{x) — ■nx/2 J;+i/2(a;). 
The Bessel transform of an oscillator state is again an oscillator state 

(^„X^) - {-iTb^nAkb''), (11) 

with a classical turning point Kn^i — Rnjb^ ~ \/4n + 21 + 3/6. This relation that can easily be derived using 
e.g. formula (7.421.4) from [3S]. The oscillator states therefore form an orthonormal basis of i^([0,oo[) that 
diagonalizes the Bessel transform. Similar properties hold for the Cartesian oscillator state based on Hermite 
polynomials, see, for example, [3]. 

An other important property that will be used in the remainder of the text is the following: 

Proposition 2.1. Let ijji a function that behaves as in x ^ 0. The projection on the oscillator state (pn,i 
can be calculated with either ipi or its Bessel transform ipi as: 

poo poo 

Cn,i= / ^pn,i{x)'ipi{x)dx^{-l)''^b Lpn^i{x)%i>i{x)dx . (12) 
Jo Jo 

Proof. Since Parseval's theorem holds we can calculate c„ either with (p^^i or its Bessel transform 

/•oo i*oo 

Cn,; = / ipn.i{x)iJi{x)dx ^ / (find{k)ipi{k)dk . (13) 
Jo Jo 

But since (pn,i{k) = {—^)"bLpn,i{kb'^), the oscillator state in the latter integral is the same as in the first 

integral, after substitution of variables and up to a constant. □ 

This result will be used in the next section to derive the asymptotic formula for c„./ with large n. This 
symmetry between ipii^) ^J^d its Bessel transform 'il'i{k) should also hold in the asymptotic formula. 

The kinetic energy operator ^ is a tri-diagonal matrix. For the oscillator basis the non-zero elements are 

r (2z + / + §)| forj = i, 

T;'.=/°°^u(.)(-^^ + ^)v',,(.)dx= for,=.-l, (14) 

[-^{t + l){t + l+'l)^ forj=z + l. 
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For the remainder of the text we will consider only the case of zero angular momentum and drop I from the 
notation. All results, however, are valid for arbitrary I. 



2.3. The asymptotic formula 

In [3 6) an asymptotic formula was derived for the projection of a smooth radial function on the state (p„. 
The derivation uses stationary phase arguments to exploit the increase of oscillations as n — > oo. We repeat 
here the derivation of the results. 

There are two main contributions to the value of the integral 

ipn{x)ip{x)dx « /o + /fl„. (15) 

There is a first contribution, denoted Iq, of the left integration boundary near zero. There is a second 
contribution, denoted /_r„, from the point of stationary phase, which is the classical turning point where 
the oscillations stop (see Figure [T]). Here we provide only general steps of calculating this integral in the 
asymptotic region. More details can be found in |36j . 

The contribution from the classical turning point can be calculated by approximating the oscillator state ipn 
near the classical turning point by an Airy function as 



^nix) 



1/6 



Ai 



1/3 



{X - Rn) 



if |a;-i?„|<l. 



(16) 



The integral representation of this Airy function and the stationary phase approximation leads to the con- 
tribution of the turning point to the integral 



bJ^^{R^). 



(17) 



The contribution of the left integration boundary, /q, can be derived by approximating the oscillator state 
near the origin by a Riccati-Bessel function 



^Pn{x) ~ (-1) 



b V TrKn 



jo{KnX), 



(18) 



where Kn is the classical turning point of the oscillator state in the momentum space. Then the contribution 
from the origin becomes 

lo^i-irUf^mn). (19) 



And the resulting asymptotic approximation of the oscillator coefficients becomes 

1 r^7,,, . , 



bJ^^{R„ 



if n > 1. (20) 
This relation has a contribution from the turning points in the coordinate space, i?„ and the Fourier space. 



Kn. Note that it satisfies the symmetry observed in 2.1 above. 

Note that Eq. ( 20 1 as well as the similar equation in [35] does not provide the order of the approximation. 
This is one of the shortcomings that will be addressed in the current paper. 



2.4-. Scattering calculations in the oscillator representation 

We now discuss how a finite linear system can be obtained that solves for the wave function in the interaction 
region and the phase shift, describing the solution in the asymptotic region. The presentation here is based 
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on the asymptotic formula and differs from how the method was derived historically. For more detail about 
the formulation of the original J-niatrix method we refer to [9] and |14j . 

Let ipi^) be a smooth two-body radial scattering state with a bounded energy, depending on one spatial 
coordinate (can be always reduced using center-of-mass relative coordinates and spherical coordinates). Since 
it is a scattering state, the function does not go to zero as a; — ?► oo. However, its Bessel transform ip{k) goes 
to zero as k ^ oo. This means that, as n grows, the only contribution to the expansion coefficient c„ comes 
from the classical turning point in coordinate space. So, for a scattering state ip with total energy E = /2 
the solution is asymptotically ji{kr) + tan((5;)n;(fcr) where 5i is the phase shift [57] holds that 

Cn « bJ J-^(i?„) = bJ J- (UkRn) + taniSi)hi{kRn)) , (21) 
V it„ y Kn \ / 

where ji and "hi are the regular and irregular solutions of the free-particle equation. The c„ becomes the 
representation of ip on the grid of classical turning points i?„. 

The aim of a one-dimensional radial scattering calculation is to find a numerical approximation to tan((5;) 
for a given potential V. This can be achieved by writing the solution in the oscillator representation as 

J c% + jnj + tan{Si) Tin J n< N 
Cn,l = < (22) 
[ Jn.i + tan((5i) Unj n> N, 

where jn.i is the regular solution of the homogeneous three-term recurrence relation 

,n — ljn—l,l + {Tn.n + 7n,n+ljri-|-l,; = 0, (23) 

where as n >• 1 holds that j„^i = b^2/RnjiikRn). The n„ / is the irregular solution of the recurrence 
relation that goes as n„ ; = 6-\/2/i?„ Ti;(fci?„) when n becomes large. 



With the form ( 22 ) we reduce the infinite linear system to a finite dimensional problem with a set of A?^ -f 1 
unknowns {c^ i^tanSi}. Once solved we simultaneously obtain the wave function of the system, {c^;}, and 
the scattering information, tan(i5/). 

A detailed description of the construction of this linear system, known as the J- matrix, can be found in |14j . 
2.5. The hybrid J -matrix and ECS method 



In |33j the asymptotic formula ( 20 1 was used to introduce a hybrid oscillator and grid representation for the 



Schrodinger equation that is useful for scattering calculations where an asymptotic form such as (21 ) is not 
explicitly known. Such a situation appears in breakup reactions of three or more particles. 

Let tpi^) be again a smooth one-dimensional radial scattering state with a bounded energy such that 

Cn ~ byl^^iRn). (24) 

V -K„ 

holds. The c„ becomes the representation of ip on the grid of classical turning points i?„. The grid distance 
between these points becomes smaller when n increases since 

h„^Rn-Rn-i-b^^. (25) 

However, the asymptotic ^ does not necessarily need to be represented on the grid of turning points 
Another option is, for example, a regular grid of equally spaced points. 

The hybrid JM-ECS method represents the one-dimensional radial wave function as a vector U' in C^+^, 
where 

^ = {co,Ci, . . . ,CN-l,i'{RN),i'{RN + h), . . . ,i^{RN + {K~l)h)). (26) 
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The first N elements represent the wave function in the oscillator representation. While the remaining K 
elements represent ip on an equidistant grid that starts at Rn, the A^-th classical turning point, and runs 
up to Rn + {K — l)h with a grid distance h equal to the difference of the last two turning points of the 
oscillator representation h ~ Rjy — Rpj^i- It is assumed that the matching point that connects the oscillator 



to finite-difference representation corresponds to a large index N such that the asymptotic formula, (24), 
applies. 

Again, the kinetic energy operator in this hybrid representation is tridiagonal since in both finite-difference 
and oscillator representations it is tridiagonal. One should only be careful in matching both representations. 
To obtain the kinetic energy in the last point of the oscillator representation, the tridiagonal kinetic energy 



formula (14) is used. It involves a recurrence relation connecting the three terms cn-2, cn-i and cn- 



The latter, the coefficient cn, is unknown. Only iI>{Rn) is available on the grid. Using the asymptotic 



relation (24), however, it is possible to calculate the required matrix element as follows: 



{Tc)n-1 = Tn-1,N-2CN-2 + Tn-1,N-1CN-1 + TN-l,NbV^/RNMRN)- 

To calculate the kinetic energy in the first point of the finite difference grid, the second derivative of the wave 
function has to be known. To approximate the latter with a finite difference formula, one needs the wave 



function in the grid points Rn-i, Rn and Rn + h. Again it is possible to apply (24) to obtain ^{Rn-i) in 
terms of cn-i'- 

CN-i/{b^2/RN-i) - 2^P{Rn) + ^(Rn + h) 



r{RN) 



/l2 



The coupling between both representations around the matching point is sketched in Figure [2j together with 
the terms involved to determine the correct matching. 



Oscillator 


Finite Differences 


Tciv-i 




Cn-2 Cjv-i 

e e 


^{Rn) ^{Rn + h) 
• • 




V ' 




iA"(fl]v) 



Figure 2: Illustration on the calculation of the kinetic energy matrix elements that are calculated in the last point -Rjv— i of the 
oscillator representation and in the first point ijjy of the finite difference representation. To calculate T applied to a solution 
vector we need to translate the oscillator representation to the grid. Vice versa for the application of the finite difference stencil 
approximation to the second derivative. 

However, in the next section we will show that the asymptotic matching condition that is used in |33] to couple 
both representations, is only accurate to ©(i?^^), while the outer grid is accurate to order 0{h^) = 0{R~^). 
Therefore the largest error is made in the matching condition that couples both representations. 

Note that introducing an absorbing layer is easy once the oscillator representation is coupled to a grid 
representation. For example ECS is implemented by extending the grid with a complex scaled part j24j . 



3. Higher order asymptotic formula 

In this section we derive a higher-order asymptotic formula. It not only takes into account the function 
values in the turning point i?„ but also the third derivative in this point. A similar asymptotic formula was 
derived by S. Igashov in [34 , however, it did not include the contributions from the Fourier space. 
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Proposition 3.1. Let '4'{x) he a regular scattering state that behaves as x' in x — 0, that is infinitely 
differentiable and (fnix) is a reduced oscillator state, solution of Then the projection of the scattering 

state on the oscillator state can be approximated as 

Cn = ^n{x)ijix)dx ^ — U(i?„) + j + O (i?,T''/') • (27) 

This expresses the expansion coefficient in terms of the wave function and its derivatives in the classical 
turning point 



Proof. The oscillator Eq. ([9| can be rewritten as 



d^ , 2R„{x-Rn) , [x-Rn) 



dx^ 



&4 



64 



iPnix) = 0. 



(28) 



For i?„ ^ |a; — i?„| (which means either n ^ 1 or |a; — i?„| ~ 0), we can neglect the quadratic term in this 
equation and get the Airy equation with the solution 



4 ^ 1/6 



b \2Rn 



Ai 



f2Rr. 



V b^ 



1/3 



{X - Rn) 



(29) 



where the normalization is chosen such that it coincides with the oscillator state near the classical turning 
point Rn- 

This can be further improved by writing tpnix) = (1 + w„(a;)) cpn"^ and inserting this in the oscillator equation. 
The equation for Un{x) becomes 



1 " (a) I (oV 1 (0)" 



Using that 



and 



(0)' 

n 

(0)" 



V.-' ^ {2Rn/b^f'^ A\' {2Rn/by\x-Rn) 



(x - R„ 



,4\l/3 



(0) 



■;^(x-i?„)Vi°'. (30) 



(2i?„/64)'/' Ai" \(2Rnlb^Y'^ {x-Rr, 



the equation (|30|) becomes 
1 



" (0) / ( '^^r: 

2/3 



1/3 



2Rn 



Ai" 



Ai' 



2Rn 



2i?„ 



1/3 



(a; - Rn) 



1/3 



{x — Rn 



+ (^Rn{x - Rn) + Rn)^^ V^i"^ ^ " ^ ~ ^n)Vi"^- 

(31) 



In the limit n — > oo, the term with Rn{x — Rn) dominates compared to (1/2) (x — Rn)"^ and the other 
terms since both the first and second derivative of Ai(a;) around zero are bounded. We find that Un{x) ~ 
— {x — Rn) / (2Rn) is a solution of the remaining equation. 

At this time we will not consider the contributions to the integral for the boundary at zero. This contribution 
will be equal to the contribution from the Fourier space as in Eq. ( 20 1 and for scattering states and large n 
this contribution is negligible. So we can lower the integration boundary and write 

2,/ar 



ip{x + Rn)Ai{x/an) 1 - 



2Rn 



dx. 



(32) 



where a„ := (6^/2_R„)^/'^ and we have substituted the integration variable x ^ x + Rn- 
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The value of this integral will be determined by the behavior around i?„, the point of stationary phase. Since 
the function ip is infinitely differentiable we can Taylor expand 



' m=0 



2Rn 



Ai(x/a„)dx 



(33) 



All integrals in these series can be calculated explicitly using specific properties of Airy function (see page 
52 of [3HI)- 



Ai(x)x^''dx = and 
3'^k! 



Ai{x)x-^''+'^dx 



Ai{x)x-^''+^dx = 0. 



(34) 



We finally get 



^£||(*""w.)-4*'»"'(«j 



- oo ^ 



(35) 



In lowest order (in terms of l/i?„) this gives us exactly the initial relation (24 1. Including the next order 
correction we get a relation (27). For testing purposes we also include terms of the next order l/R^ 



c„ = — ( 1p{Rn 



(36) 

□ 



Note that a similar result is readily obtained for Cartesian harmonic oscillator states that are based on 
the Hermite polynomials. There, however, there is a contribution from both turning points, one from 
i?„ = y^2n + 1 and i?„ = ~y^2n+l. 

Corollary 3.2. Let ^£ function that behaves as x^ in x — that is infinitely differentiable, then the 

asymptotic expansion coefficient is 

V?„(x) 2P{x)dx = bJ--{ 7A(i?„) + ^^"'(-R„) 

V V 



(37) 



Proof. The asymptotic formula should have the same result when we interchange "0(2;) with its Bessel trans- 



form ipik) as a result of Proposition 2.1 Indeed, the integral 



c„ = (-1)"6 / ^„{k)^{k)dk 
Jo 

will have a contribution from integral boundary near and a contribution from the turning point Kn in 
Fourier space. The latter is, with the help of the previous rcsuhs, {-lYby/2jKn{'ip{Kn) + 'i}^'" (Kn) /Qb'^Kn). 
While the contribution from boundary become the contribution from the turning point i?„ in coordinate 
space. 

□ 



Example 3.3. We illustrate the convergence of the asymptotic formula (371 with an application to the 
function ip{x) — exp(— ax) sin(a;). We give results for two choices of a. As the scale of the representation 
is defined by the oscillator length b, the result will change with its choice. If the product ab is large, the 
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n 



Figure 3: Convergence of the general asymptotic formula ( |37| with only coordinate terms (crosses), only Fourier terms (circles), 
both coordinate and Fourier terms (solid line) for test function t/'(x) = exp(— ax) sin(x) with different values of a: left figure — 
a = 0.2, right figure — a = 1.5 (b = 1 in both cases) 



Fourier components will dominate in the expansion coefficient. When ab is small the function ip resembles 
a scattering state, on a scale defined by b, and the coordinate approximation will dominate. However, as 
Figure^ illustrates the combined formula always gives the correct result. We see an overall convergence with 
r/iis is smaller than Oiji^) which goes as 0{l/n). In general, it is possible to construct a function, 
for which both coordinate and Fourier terms fail to represent the exact oscillator coefficient, but the combined 
formula remains valid even in this case (see Figure^^. 

Example 3.4. We have verified the asymptotic relations derived in the previous section with a numerical 
experiment using the wave function ip{x) = sin(fca::). We calculate oscillator coefficients directly and compare 
them against the values obtained with asymptotic formulas of different order derived previously (see Fig.lSV. 



3.1. Inverse relation 

The asymptotic expansion coefficient is an approximation with the help of the functions values of tp evaluated 
at certain grid points. It is also useful to derive in inverse relation that allows us to construct the function 
value in certain grid points given the expansion coefficients. 



To approximate c„ with Eq. (271 the function value and its third derivative need to be calculated at the 
turning point When only the function values of ip are available at the turning points i?„, the third 
derivative can be approximated by finite differences. 

The coefficient is then calculated as 

cn = HRn) + ^ E Df^mn+k)^ + 0{R-''l^) + K^'^e^, (38) 
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Figure 4: Convergence of the general asymptotic formula jSTjl for the test function tp{x) = exp(— ax) cos(x) with a = 0.2, for 
which neither coordinate (crosses) nor Fourier terms (circles) give accurate result, while the combined formula (solid line) does. 




Figure 5: Convergence of the asymptotic formula (jSGjI (left) and the inverse asymptotic formula l |41| (right) with different 
number of terms included (crosses — one term, solid line — two terms, circles — three terms). The test function has the form 
■(/>(a;) = sin(fex) with k = 1.5 and the oscillator length b = 0.8 
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where D^^^^ is the matrix with coefficients that approximates the third derivative and k indicates the stencil 
points and e„ is the error term of this approximation. In case of 5-point finite difference approximation 
k e {—2, —1, 0, 1, 2}. Note that the grid of turning points i?„ is an irregular grid and the coefficients will 
depend on the local distances between the neighboring grid points. The error of the approximation e„ will 
also depend on the local grid distances as 0{hP) ~ 0{R~p), where p is the order of the approximation. 

The resulting transformation can be represented by a banded sparse matrix U 

Cn^Y. Unk^iRk) + 0{R-'''^). (39) 

k 

Note that the matrix elements of U will only depend on the values of i?„. Indeed, The differential operator 
D^'^) only depends on the grid distances and so do the coefficients of the asymptotic expression. 



The relation (38 1 translates if) on the grid of turning points to corresponding c„. It is also useful to derive the 



inverse relation that gets ipiRn) from known values of c„. We can not use a direct inversion of Eq. (27) since 



it involves the third derivative, a dense operator in the oscillator representation. However, an approximate 



inverse relation can be obtained by rearranging terms in ( 38 1 as 



"^(^"^ = y T^^" - ^ E D^k ^^(Rn+k) + 0{R-') + R-hn (40) 
and replacing values of the wave function in the right-hand side by oscillator coefficients using only the first 



term of (24|, i.e. ipiRn) = (1/6) v/^„/2 c„ + 0{R-^). This only introduces an error of the order of 0{Rj^ ^) 



but combined with the l/i?„ this gives the approximate inverse relation 



fc 

that is accurate to 0{R^^). An example of the resulting numerical accuracy of this inverse asymptotic 
relation is shown on the right panel of figure |5] Also this transformation can be presented as a sparse banded 
matrix multiplication 

V'(i?„) = ^^^"^ + ^(^n')- (42) 

n 

Again the matrix elements of W only depend on the values of the turning points i?„ . 
Combining the two relations leads to an approximate partition of unity: 

UW ^I + 0[R-^). (43) 

It is important to note that since these transformation matrices only depend on values of i?„ they can, in 
principle, both be defined for an arbitrary grid. 



3.2. Approximate discretization of Operators 

With the help of these two transformation matrices U and W we can now build an approximate oscillator 
representation of an operator Q using its finite difference representation. Let Q^^'^'> be the finite difference 
representation of Q on the grid R„. Then the application of the Q on -0 can be written as 

iQi^)ir^ = E = E + ^(""')' ^here Qt^^ := [UQ^^^'^WU^. 

m m 

We illustrate the accuracy of this relation on a second derivative operator D'-2), as we intend to apply the 
constructed representation to Hclmholtz-type equations. 
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Figure 6: Convergence of the approximate second derivative operator D = UD^-f'^^W (crosses) and the approximate identity 
operator / = UW (circles) with a test function of the form ^{x) = sin(fca;) (k and b are the same as in Figure|5]l 



£)(2)(osc) ^ [/£)(2)(/rf)|y-^ (44) 

where jg g, finite difference matrix of the second derivative. To analyze the accuracy of this approx- 

imation Figure [g] shows the result of the error operator defined as (i5(°*'=) — acting on the vector of 

oscillator coefficients of the test function 'ip{x) = sin(fcx). For comparison we show the error of the approx- 
imate identity operator, that can be defined as (/ — /) — (/ — UW). We see that both considered error 
operators decay as 0{n^^) in high-n region. This means that our approximate operators are asymptotically 
equal to the exact ones and we can expect the same order of convergence in a solution of the scattering 
problem. 



It is important to note that all matrices in ( 44 1 can be built for any arbitrary spatial grid, though it will 



not be an approximate oscillator representation. But we can modify this grid only in the asymptotic region 
{e.g. to implement the absorbing boundary with ECS). Then the coefficients c„ have no physical meaning 



as oscillator coefficients but the coordinate wave function can sill be reconstructed with ( 42 1 



4. High-order hybrid representation for scattering problems 

The aim is now to build a hybrid representation for scattering calculations where the oscillator basis is used in 
the internal region and finite differences in the outer region. The matching of the two should use the higher- 
order asymptotic formula (27). However, the strategy displayed on Fig. [2] cannot be easily applied since the 
third derivative of the wave function needs to be estimated from several neighboring points symmetrically 
distributed on both sides of the matching point. But in the matching point we only have the required 
function values on one side of the matching point. This arguments holds both for the oscillator and for the 
finite difference representation. We therefore present a matching strategy that differs from the proposal of 



Consider a arbitrary grid i?„ in [0, oo] with n = 1, 2, . . . that differs from the grid of turning points i?„. On 
this grid we can construct the differential operator D^'^' and construct the operators U and W by replacing 
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every appearance of i?„ in equations (38) and (41 1 by i?„. Similarly, for a function tjj sampled in the grid 
points Rn we can calculate 



k 



Unki'iRk)- 



(45) 



These coefficients c„ are not the expansion coefficients of the function %p in the oscillator function basis. Only 
when Rn equals i?„, the oscillator turning points, then the c„ are approximations to the oscillator expansion 
coefficients c„. 

To build the hybrid representation, we choose the grid Rn such that first + 1 points correspond to the 
turning points The remaining points of i?„, for n > + 1, are chosen to correspond to points on a 
equidistant finite difference grid. To solve a scattering problem the equation ^ needs to be discretized with 
finite differences on the grid Rn and then transformed with the help of tj and W into 



E 

i 



(46) 



to arrive at an equation for c,;. 

The first N coefficients of c„ correspond now to approximations to oscillator expansion coefficients c„. 
However, because n is low, they are only a poor approximation. The idea is now to replace the first N 
coefficients with the exact coefficients. At the same time we replace the first N rows of the matrix with the 
exact operator in the oscillator representation. The linear system then becomes 



„(osc) 
-'^10 


u-(osc) 

^OiV 

rr(osc) 


tt{0Sc) 

-"OJV + 1 
„(osc) 
-"lJV + 1 


rr(osc) 


-"iVJV ^ 


tt(oSc) 

-"nn+1 


[i7//(^'')iy]jv+i,o . 




[U{H(f''^ - E)W]m+i,n+i ■•• 



J 



Co 
ci 



C]v 



cjv+i 



V 



J 



(osc) 


(osc) 



p{osc) 

■In 



V 



(47) 

where H'^"^'^^ is the representation of the Hamiltonian in the oscillator representation and H'^^'^^ in the finite 
difference representation. 

We emphasize the difference with [33; . Here we do not match two regions by using the asymptotic formula in 
one point. Now the representation in the asymptotic region is a fairly good approximation of the oscillator 



representation. Therefore the structure of the discretized wave function is simpler than ( 26 ) . Now we have 

(cq, . . . Cjv, Cjv+1, . ■ • ck) , 



where in the initial version of hybrid method, the latter elements of the vector are the function values of ip in 
the grid points. Now the internal region is covered by an exact oscillator representation, and the asymptotic 
region is covered by approximate representation which is based on finite differences and includes the ECS 
transformation. 



Numerical Illustration 



We illustrate the method for a one-dimensional radial example. After the solution of (47 1, we first reconstruct 
the coordinate wave function using (42 1. Outside the range of the potential V the solution can be written as 
a linear combination ipse = Ahf {kx) + BK[ [kx) , where li^ are the in- and outgoing Riccati-Hankel functions. 
The coefficient A is then extracted with 



(V',c(x), V(M) /W {h+{kx),h-{kx)) , (48) 



where x is outside the range of the potential but still on the real part of the ECS domain. The Wronskian 
is calculated as W{u, v) = u'v — v'u, where the derivatives can be implemented with finite differences. From 
the Wronskians for A and B we can extract the phase shift of the solution. 
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Figure 7: Scattering pliase shift of tlie model problem of Section 4 with a Gaussian potential (left) and the absolute error of 
the scattering phase shift calculated with JM-ECS (dashed line) and our new approach (solid line) depending on the energy of 
the system (right). Both calculations were made with 100 functions in the oscillator basis with the oscillator length b = 0.7. 




Figure 8: Accuracy of the scattering phase shift calculated with JM-ECS (crosses) and the new method (circles) depending on 
the size of the oscillator basis. Calculations were made for two values of energy: E = 0.2, where S = 43.73° (left) and E = 3, 
where S = 19.67° (right) 



It is important to note that the accuracy of the method does not directly depend on the size of the asymptotic 
region as long as this region is located outside the range of the potential {V ^ in the considered part of 
H'^"'"^). In this region we use a grid of 150 equidistant points and then the ECS layer that spans 10 
dimensionless length units and applies a complex rotation of 45° to the coordinate axis. 

We first consider a model problem with a attractive potential in Gaussian form, V{x) — — exp(— x^), and we 
limit our consideration to the case of zero total angular momentum I = 0. Nevertheless, all the conclusions 
are applicable to higher angular momenta as well. The right panel of Figure [7] shows the error in the 
scattering phase shift of the considered problem calculated with the original JM-ECS and the new approach. 
To obtain the reference phase shift we use the highly accurate variable phase approach (VPA). We see that 
for all energies the new method gives much more accurate and less oscillatory results. The convergence as a 
function of the number of oscillator states in the inner region of both methods is shown on Figure [8j We see 
that we arrived at the desired the convergence rate to N^^ in both low-energy and high-energy regions. 
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Figure 9: Same as Figure [s] but for the Morse potential (two upper panels) and Yukawa potential (two lower panels) 

The Gaussian shape of the interaction potential was chosen as the main model problem as it is particularly 
well adapted to the harmonic oscillator basis and J-matrix method for this potential converges for any 
oscillator length b. However we can also test our approach against other short-range potentials. But if 
the asymptotic behavior of the potential is not Gaussian then the results will converge only for specific 
values of the oscillator length, that match the range of the potential. Also, if the potential has a special 
point in the origin (like Yukawa potential) then the convergence will be much slower due to the importance 
of the Fourier contribution in the potential matrix elements. For additional convergence tests we have 
chosen two potentials with non Gaussian asymptotics and different behavior in the origin: Morse potential 
V{x) = exp(— 2|a;/6|) — exp(— and Yukawa potential V{x) — — exp(— |a;/6|)/a;. In both potentials the 
potential range is chosen to match the oscillator length of the basis. Figure |9] shows the convergence of the 
phase shift for these two potentials. We see that in high-energy region the convergence pattern is mostly 
similar for all potentials as the kinetic energy is much higher than the potential energy in this case. The 
behavior of the error is much less clear for low energies. The convergence rate appears to be the same here 
for both methods with Morse potential, for Yukawa potential the convergence rate is blurred due to high 
oscillations which clearly indicates the importance of Fourier terms. Though we can generally conclude that 
the new approach always gives better and generally less oscillatory result. 
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5. Extension to systems with more particles 



The results of the previous sections can be generahzed to systems of three and more particles. Let ri and 
r2 be two relative coordinates describing a three-body problem. The 6D wave function can be expanded in 
partial waves 

vE'(ri,r2) = *(pi,p2,ei,^2,0i,02)= V "^'"'""""^ Yi^,,^, (02, 02). (49) 

, , PlP2 

ti ,mi ,t2 ,m2 

Such an expansion is used for example to describe double ionization processes in atomic physics |39j 

Let ^^{x, y) now be an infinitely differentiable two-dimensional radial scattering wave function that is ex- 
panded in the bi-oscillator basis as 



oo oo 



tp(x,y) Cnmfn{x)(Pm{y), (50) 



n— m— 



where x and y should be interpreted as two radial coordinates. The expansion coefficient is then calculated 
as a double integral that can be approximated by applying the asymptotic relation twice. First in the 
x-direction and then in the ^/-direction 



oo roQ 



ipn{x)(prn{y)il^{x, y) dxdy 



dx 



= 2{RmRn) ^^"^i^iRn, Rrn) H Rn^ Rm d^xxi^iRn, R-m) H Rn^^^ Rm^^''Pyyy{Rn, Rm) 

(51) 

It is important to note that the accuracy depends on both indices n and m. For example, when the low order 
approximation is used for both integrals, there will be error terms of the order RVri^^R^^ and Rn^^"^ R^'^ . 
Along the diagonal, i.e. when i?^ = Rn, it is accurate up to order Rn^^'^R'^^^ = R^^ — 0{n^^), which is of 
the order h^, with /i„ the distance between the classical turning points. For the higher-order approximation 
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the error term along the diagonal, where i?„ ~ will be R^^ = 0{n^^/'^). In the left panel of Figure 
we compare the exact expansion coefficients with the first order and second order approximation along the 
diagonal for n = m. The results show an improvement with the higher-order expression over the first order 
approximation . 

However, when one of the indices remains small, for example m, then the error becomes asymptotically 

— 1/2 i 1 

0{Rn ), as expected from (51 1. As n becomes large, all error terms with a higher-order accuracy decay 
until the term 0{Rn^^'^ K^^^) dominates for the higher-order formula, or 0{Rn^^^ R„?^^) for the low order 
accuracy. This is shown in the right panel of Figure |10[ 



5.1. Three-body scattering problems 

In a similar way as in the previous sections we can use the asymptotic formulate to build a hybrid repre- 
sentation that can solve the radial equation of two radial variables that arises after an expansion of a full 
three-body scattering problem equation in spherical harmonics ( 49 ) . This typically leads to a set of coupled 
two-dimensional radial partial waves. A diagonal block of this coupled set reads 

Id' I d^ h{h + l) hih + l 
' 2 dx2 2 dy^ 2a;2 + 2y^ 



+ V{x,y)-E]^P{x,y)^x{x,y) (52) 
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n,m n 

Figure 10: Convergence of the asymptotic formula for two-dimensional radial functions. The error is calculated by comparing 
the exact expansion coefficient with the asymptotic approximations. The function f{x, y) = sin(fca::) s'm{ky) along the diagonal 
(n=m, left) and along the line with fixed m=20 (right). Crosses show the first order and solid lines correspond to the second 
order asymptotic expansion coefficient (fc = 1, = 1). 



with boundary conditions ip{x, 0) = ''PiO, y) ~ 0. 



The solution of this equation can be represented in a bi-oscihator basis, (50), using two indices n and m. 
Then similar to the two-body case we use approximate oscillator coefficients in the asymptotic region. But 
there are now three asymptotic regions: first, the region where n is large, second, the region where m is 
large and finally, the region where both n and m are large. For each of these regions we define approximate 
oscillator coefficients that take the form: 



Cm„:=Vt/„W ip„^iy)i^iR^,y)dy + OiRJ^^) where n>l (53) 
Jo 

Cnm:=Vt/mi/ (pn{x)tp{x, R,) dx + 0{R;^^/^) whcrc m>l (54) 
Jo 

Cnm ■■= J2 UmUrnMR^^ ^j) + O {R-^ ' ^ R-^ / ^) + O {R-^^ ' ^ R-'> ' ^) whcrc u, m » 1, (55) 



where the error terms in the last expression are explained in (51 ). 

Then the coefficient matrix Cnm of hybrid representation is divided into four blocks corresponding to different 
regions 



Coo 




Com 


■ CoK' 


ClO 


■ CiM-1 


C\ M 


■ ClK' 


Cat-IO • 


■ Cjv_iM-l 


cn-im ■ 


■ Cn-IK' 


cjvo 


■ Cn M~i 


Cn m 


■ Cn K' 


cjv+10 • 


■ Cn+IM-1 


cn+im ■ 


■ cn+ik' 


cko 


■ CkM-1 


Ck M 


■ Ck K' 



(56) 



/ 

where we define N and M as sizes of the oscillator bases in each direction and K and K' are the total 
number of variables in each direction. 
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To build a linear system corresponding to ( 52 ), we reshape the matrix c„m to a vector. Then the Hamiltonian 



matrix of a 2D problem is constructed as a Kronecker sum of the two-body Hamiltonian's constructed as 



in ( 47 ) and a two-body potential matrix. In the case of our approximate oscillator representation this takes 



the form 

i, 



so the Kronecker sum contains the approximated unity operator instead of the real one. Where and 
are the transformation matrices for the x coordinate and similarly for the y coordinate. The total size of the 
Hamiltonian matrix is K x K' . 

5.2. Numerical results 

We first evaluate the method with the help of a model Helmholtz problem with a constant wave number or 
energy E = fc^/2 and for li = and I2 = 0. The equation is then 

-^A~E^tP{x,y)^^o{x)My), (58) 

with boundary conditions i'ix, 0) = and ■0(0, u) — 0. The form of the right-hand side here was chosen 
for simplicity as it results in a vector (1, 0, 0, ... 0) in bi-oscillator representation. This problem is exactly 
solvable with the help of the Greens function 

Gix,y;x',y') =- [4'\Wix - x'Y + {y - y'f) - H^^\k^{x + x'f + [y - y'f) (59) 
- Hl^\k^{x-x'Y + {y + y'Y) + H^^\k^ {x + x'Y + {y + y'Y)) , (60) 

where Hq^^ is a 0-th order Hankel function of the first kind. The scattering solution is then 



00 /'OO 

I l\ I I J I J I 



i'{x,y)^ \ \ LpQ{x')ifQ{y')G{x,y;x' ,y') x'y' dx' dy' (61) 
Jo Jo 

For simplicity we are not going to compare any scattering information extracted from the wave function as 
it involves additional operations like surface integration, which can lead to additional loss of accuracy and 
need to be studied separately. We will compare the values of the wave function in a fixed spatial point. As 
the spatial grid is different for every size of the oscillator basis we can not ensure that any fixed spatial point 
lies exactly on the grid point at every calculation, so we need to interpolate. We have used a cubic spline 
interpolation, and from comparison to the other interpolation methods we can expect that the additional 
error introduced by this operation is negligible compared to the errors of the method (of course, that is only 
if the function is smooth enough, that means not very high values of E). 

As we do not have any potential in this model problem, the convergence of the proposed hybrid oscillator 
representation will depend mainly on the accuracy of the asymptotic relations used. This means that we 
can choose the size of the spatial domain on which we solve the two-dimensional radial problem as small as 
possible but not smaller than the region spanned by the exact oscillator basis. As the maximal size of the 
oscillator basis we use in this calculations is 350 and the oscillator length was chosen as 6 = 0.7 {Rm « 26) 
the size of the spatial domain was chosen as 30 dimensionless length units of the real grid and additional 10 
units of the ECS layer. The value of the wave function was extracted at the point (28, 28) when we increase 
the basis size simultaneously and at (28, 12) when we keep the basis size fixed in the y-dimension. 



In Figure 11 we compare the numerical solution of the linear system with the one obtained from (61 ) which 



we consider as exact. We see that, similar to Figure 10 if we expand the basis only in one direction the 
convergence is of the order N~^^'^. 
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Figure 11: Convergence of the two-dimensional scattering wave function for the problem in Eq. \58\ at the fixed spatial point. 
Left figure represents the case where the basis was increased in two dimensions simultaneously {N = M), right figure shows the 
convergence with a fixed M = 50. Both calculations are made with b = 0.7 and E = 2. 



Next we test the method on a model potential scattering problem described by ( 52 ) with zero partial angular 
momenta and the Gaussian interaction potential in the form 

V{x, y) = -3e-^' - Se"^' + lOe-^'-^' (62) 

The one-body potential, here Vi{x) — — 3exp(— x^), can support one bound state. This means that using 



the potential ( 62 ) we can model a three-body breakup problem using the right-hand side in the form 

X{x,y) = -{V{x,y) - Vi{x)) h{x)Uv) . (63) 

where /o(a;) is the wave function of the bound state, jo(2/) is the Riccati-Bessel function and together they 
represent the initial state of the system. For the considered problem there are no analytical results to 
compare with, so on the Figure [12] we present the solution of the linear system for Cnm and the spatial wave 
function, reconstructed from c„m by applying the transformation matrix W ® W. On both figures we can 
clearly recognize the patterns of elastic scattering (the plane waves along the axes) and the breakup (the 
radial waves with lower frequency as part of the energy was spent to break the bound state). The in-depth 
analysis of the three-body scattering results is the subject of future work. 



6. Conclusions and outlook 

This paper focuses on scattering calculations in the oscillator representation, where the solution is expanded 
in the eigenstates of the harmonic oscillator. The oscillator representation is not the most natural represen- 
tation to describe scattering processes since it involves a basis set that is designed to describe bound states. 
This leads to a rather low convergence and may result in a linear system with very large dense matrix. It is 
often more natural to describe a scattering process with the help of a grid representation. These grid-based 
calculations have been very successful in describing scattering and breakup processes in atomic and molecular 
physics. The Helmholtz equation is also efficiently solved on these grids. 

However, internal structure of the nuclear clusters and other many-particle systems are efficiently described 
by such an oscillator basis, since the bottom of the potential can be mimicked by the oscillator potential 
leading to an efficient description of the internal structure. 
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Figure 12: The solution of tlirce-body scattering problem in oscillator representation (left figure) and the reconstructed spatial 
wave function (right figure). The calculation is made with b = 1 and E = 1. The problem has Dirichlet boundaries on all sides. 

In this paper we combine the advantages of grid-based calculations with those of the oscillator representation. 
The method was originally proposed in [33! , as a further development of the J- matrix or algebraic method for 
scattering. There the method combined the grid and oscillator representation with a low-order asymptotic 
formula. In this paper we have improved this matching with a higher-order approximation and this brings 
the overall error down to the level of the discretization error of the finite-difference grid. 

Although a similar asymptotic formula appeared earlier in the work of S. Igashov [34) . we believe that the 
asymptotic formula presented in the current paper is more generally applicable. Furthermore, we have used 
the asymptotic formulas to improve the accuracy and convergence of the hybrid simulation method. The 
convergence of the method is illustrated with various examples from two-body and three-body scattering. 

In preparation of this papers, our initial efforts involved application of the strategy showed in Figure [2] with 
the higher-order formula. However, this required the use of a forward and backward stencils to estimate the 
third derivative in the first grid point of the finite difference representation. This strategy, however, gives rise 
to eigen modes with a negative energy localized near the interface. This destroys the positive definiteness 
of the Laplace operator. These modes are avoided in the proposed method that uses a symmetric stencil 
around the matching point. This gives satisfactory convergence. 

We have applied the asymptotic technique to radial oscillator state that are based on Laguerre polynomials. 
Similar results can be derived for the ID oscillator states that are based on the Hermite polynomials. 
These Hermite polynomials are closely related to Gaussian Type Orbitals (GTO) that are frequently used in 
computational quantum chemistry |40j . It is worth to explore if the asymptotic formulas make it possible 
to combine GTO's with grid based calculations to describe molecular scattering processes. 

In the future further research is necessary on asymptotic expressions of products of two functions. This will 
involve convolution integrals. Better expressions for products will also help to take into account asymptotic 
behavior of the potentials in the grid representation. 
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